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Convection in the solar interior is thought to comprise structures on 
a spectrum of scales. This conclusion emerges from phenomeno- 
logical studies and numerical simulations, though neither covers the 
proper range of dynamical parameters of solar convection. Here, 
we analyze observations of the wavefield in the solar photosphere 
using techniques of time-distance helioseismology to image flows in 
the solar interior. We downsample and synthesize 900 billion wave- 
field observations to produce 3 billion cross-correlations, which we 
average and fit, measuring 5 million wave travel times. Using these 
travel times, we deduce the underlying flow systems and study their 
statistics to bound convective velocity magnitudes in the solar inte- 
rior, as a function of depth and spherical-harmonic degree t. Within 
the wavenumber band t < 60, Convective velocities are 20-100 times 
weaker than current theoretical estimates. This suggests the preva- 
lence of a different paradigm of turbulence from that predicted by 
existing models, prompting the question: what mechanism trans- 
ports the heat flux of a solar luminosity outwards? Advection is 
dominated by Coriolis forces for wavenumbers t < 60, with Rossby 
numbers smaller than ~ 10 -2 at t/Rq = 0.96, suggesting that 
the Sun may be a much faster rotator than previously thought, and 
that large-scale convection may be quasi-geostrophic. The fact that 
iso-rotation contours in the Sun are not co-aligned with the axis of 
rotation suggests the presence of a latitudinal entropy gradient. 

Introduction 

The thin photosphere of the Sun, where thermal transport is domi- 
nated by free- streaming radiation, shows a spectrum in which gran- 
ulation and supergranulation are most prominent. Observed proper- 
ties of granules, such as spatial scales, radiative intensity and pho- 
tospheric spectral-line formation are successfully reproduced by nu- 
merical simulations [1,2]. In contrast, convection in the interior is 
not directly observable and likely governed by aspects more diffi- 
cult to model, such as the integrity of descending plumes to diffusion 
and various instabilities [3]. Further, solar convection is governed by 
extreme parameters [4] (Prandtl number ~ 10 ~ 6 — 10 -4 , Rayleigh 
number - 10 19 - 10 24 , and Reynolds number - 10 12 - 10 16 ), which 
make fully resolved three-dimensional direct numerical simulations 
impossible for the foreseeable future. It is likewise difficult to repro- 
duce them in laboratory experiments. 

Turning to phenomenology, mixing-length theory (MLT) is pred- 
icated on the assumption that parcels of fluid of a specified spatial 
and velocity scale transport heat over one length scale (termed the 
mixing length) and are then mixed in the new environment. While 
this picture is simplistic [5], it has been successful in predicting as- 
pects of solar structure as well as the dominant scale and magnitude 
of observed surface velocities. MLT posits a spatial convective scale 
that increases with depth (while velocities reduce) and coherent large 
scales of convection, termed giant cells. Simulations of anelastic 
global convection [6, 7, 14, 15], more sophisticated than MLT, sup- 
port the classical picture of a turbulent cascade. The ASH simulations 
[6] solve the non-linear compressible Navier-Stokes equations in the 
anelastic limit, i.e., where acoustic waves, which oscillate at very dif- 
ferent timescales, are filtered out. Considerable effort has been spent 
in attempting surface [8] and interior detection [9, 10] of giant cells, 
but evidence supporting their existence has remained inconclusive. 

Results 

Here, we image the solar interior using time-distance helioseismol- 
ogy [9, 10, 11]. Raw data in this analysis are line-of- sight photo- 



spheric Doppler velocities measured by the Helioseismic and Mag- 
netic Imager [12] onboard the Solar Dynamics Observatory. Two- 
point correlations from temporal segments of length T of the ob- 
served Doppler wavefield velocities are formed and spatially aver- 
aged according to a deep-focusing geometry [13] (Figures and ). We 
base the choice of T on estimates of convective coherence timescales 
[16, 17, 6]. These correlations are then fitted to a reference Gabor 
wavelet function [18] to obtain travel-time shifts 5t(0, 0, T), where 




Fig. 1. Line-of-sight Doppler velocities are measured every 45 seconds at 
4096 X 4096 pixels on the solar photosphere by the Helioseismic and Magnetic 
Imager (background image). We cross correlate wavefield records of temporal length 
T at points on opposing quadrants (blue with blue or red with red). These "blue" 
and "red" correlations are separately averaged, respectively sensitive to longitudinal 
and latitudinal flow at (6,cf);r/R@ = 0.96), where (0, <p) is the central point 
marked by a cross (see Figure for further illustration). The longitudinal measurement 
is sensitive to flows in that direction while the latitudinal measurement to flows along 
latitude. We create a travel-time maps Sr(6, (f),T) by making this measurement 
about various central points (9, cf)) on the surface. Each travel time is obtained upon 
correlating the wavefield between 600 pairs of points distributed in azimuth. 
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(0, (/>) are co-latitude and longitude on the observed solar disk. By 
construction, these time shifts are sensitive to different components 
of 3D vector flows, i.e., longitudinal, latitudinal or radial, at specific 
depths of the solar interior (r / Rq = 0.92, 0.96) and consequently, 
we denote individual flow components (longitudinal or latitudinal) 
by scalars. Each point (6, <j)) on the travel-time map is constructed 
by correlating 600 pairs of points on opposing quadrants. A sample 
travel-time map is shown in Figure . 

Waves are stochastically excited in the Sun, because of which the 
above correlation and travel-time measurements include components 
of incoherent wave noise, whose variance [19] diminishes as T _1 . 
The variance of time shifts induced by convective structures that re- 
tain their coherence over timescale T does not diminish as T _1 , al- 
lowing us to distinguish them from noise. We may therefore describe 
the total travel-time variance cr 2 (T) = ^2 e ^{Sr 2 (0, 0, T)) as the 

sum of variances of signal S 2 and noise N 2 /T, assuming that S and 
N are statistically independent. Angled brackets denote ensemble 
averaging over measurements of 5r(6, 0, T) from many independent 
segments of temporal length T. Given a coherence time T co h, we fit 
cr 2 (T) = S 2 + N 2 /T over T < T coh to obtain the integral upper 
limit S. The fraction of the observed travel-time variance that can- 
not be modeled as uncorrected noise is therefore S 2 / a 2 (T co h.) . For 
averaging lengths T co h (= 24 and 96 hours) considered here, we find 
this signal to be small, i.e., S 2 <C Af 2 /T co h, which leads us to con- 
clude that large-scale convective flows are weak in magnitude. Fur- 
ther, since surface supergranulation contributes to S, our estimates 
form an upper bound on ordered convective motions. 

Spatial scales on spherical surfaces are well characterized in 
spherical harmonic space: 

Sr em (T)= f sinOdO f dtp 5t(Q, 0, T) Yl m (6, 0), [1] 
Jo Jo 

where Yim are spherical harmonics, (£, m) are spherical harmonic 
degree and order, respectively, and Sre m (T) are spherical harmonic 
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Fig. 2. The cross-correlation measurement geometry (upper panel; arrowheads - 
horizontal: longitude, and vertical: latitude) used to image the layer r / Rq = 0.96 
(dot-dashed line). Doppler velocities of temporal length T measured at the solar sur- 
face are cross correlated between point pairs at opposite ends of annular discs (colored 
red and blue); e.g., points on the innermost blue sector on the left are correlated with 
diagonally opposite points on the outermost blue sector on the right. Six-hundred cor- 
relations are prepared and averaged for each travel-time measurement. Travel times 
of waves that propagate along paths in the direction of the horizontal and vertical 
arrows are primarily sensitive to longitudinal and latitudinal flows, and vq, re- 
spectively. The focus point of these waves is at r / Rq = 0.96 (lower panel) and 
the measured travel-time shift Sr(6, (j), T) is linearly related to the flow component 
v(r I Rq = 0.96, </>) with a contribution from the incoherent wave noise. We 
are thus able to map the flow field at specific depths v(r, 6, (j>) through appropriate 
measurements of Sr(0 : (f), T). For the inversions here, we create travel-time maps 
of size 128 X 128 (see Figure ). For reference, we note that the base of the con- 
vection zone is located at r / Rq = 0.71 and the near-surface shear layer extends 
from r / Rq = 0.9 upwards. 



coefficients. Here, we specifically define the term "scale" to denote 
2ttRq I \J £{£ + 1), which implies that small scales correspond to 
large £ and vice versa. Note that a spatial ensemble of small con- 
vective structures such as a granules or inter- granular lanes (e.g., as 
observed on the solar photosphere) can lead to a broad power spec- 
trum that has both small scales and large scales. The power spec- 
trum of an ensemble of small structures, such as granulation patterns 
seen at the photosphere, leads to a broad distribution in £, which we 
term here as scales. Travel-time shifts SStm, induced by a convec- 
tive flow component vim(r), are given in the single-scattering limit 
by SSim = J r 2 dr JCi (r) V£ m (r), where Kg is the sensitivity of the 
measurement to that flow component. The variance of flow-induced 
time shifts at every scale is bounded by the variance of the signal in 
observed travel times, i.e., (5S 2 &m ) < S 2 /a 2 (T coh ) (Sr^ m (T coh )) . 
To complete the analysis, we derive sensitivity kernels /Q(r) that al- 
low us to deduce flow components in the interior, given the associated 
travel- time shifts (i.e., the inverse problem). 

The time-distance deep-focusing measurement [13] is calibrated 
by linearly simulating waves propagating through spatially small flow 
perturbations, implanted at 500 randomly distributed (known) loca- 
tions, on a spherical shell at a given interior depth (Figure ). This 
delta-populated flow system contains a full spectrum, i.e., its power 
extends from small to large spherical harmonic degrees. The simu- 
lated data are then filtered both spatially and temporally in order to 
isolate waves that propagate to the specific depth of interest (termed 
phase- speed filtering). Travel times of these waves are then measured 
for focus depths the same as the depths of the features, and subse- 
quently corrected for stochastic excitation noise [22]. Note that these 
corrections may only be applied to simulated data - this is because 
we have full knowledge of the realization of sources that we put in. 
Longitudinal and radial flow perturbations are analyzed through sep- 
arate simulations, giving us access to the full vector sensitivity of this 
measurement to flows. Travel-time maps from the simulations appear 
as a low-resolution version of the input perturbation map because of 
diffraction associated with finite wavelengths of acoustic waves ex- 
cited in the Sun and in the simulations. The connection between the 
two maps is primarily a function of spherical-harmonic degree £. To 
quantify the connection, both images are transformed and a linear re- 
gression is performed between coefficients of the two transforms at 
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Fig. 3. A travel-time map consisting 16,384 travel-time measurements, spanning a 
60° X 60° region (at a resolution of 0.46875 degrees per pixel) around the solar disk 
center, obtained by analyzing one day's worth of data taken by the Helioseismic and 
Magnetic Imager instrument [12] onboard the Solar Dynamics Observatory satellite. 
3.2 billion wavefield measurements were analyzed to generate 10 million correlations, 
which were averaged and fitted to generate this travel time map. This geometry and 
these particular wave times are so chosen as to be sensitive to flow systems in the 
solar interior. The spectrum of these travel times shows no interesting or anomalous 
peaks that meet the detection criteria (described subsequently). 
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each £ separately (see online supplementary material for details). The 
slope of this linear regression is the calibration factor for degree £. 

We apply similar analyses to 27 days of data (one solar rota- 
tion) taken by the Helioseismic and Magnetic Imager from June- July 
2010. These images are tracked at the Carrington rotation rate, inter- 
polated onto a fine latitude-longitude grid, smoothed with a Gaussian, 
and resampled at the same resolution as the simulations (0.46875 
deg/pixel). The data are transformed to spherical harmonic space 
and temporal Fourier domain, phase- speed filtered (as described ear- 
lier) and transformed back to the real domain. Cross correlations and 
travel times are computed with the same programs as used on the 
simulations. Strips of 13 deg of longitude and the full latitude range 
are extracted from each of the 27 days' results and combined into 
a synoptic map covering a solar rotation. The coefficients from the 
spherical harmonic transform of this map are converted, at each de- 
gree £, by the calibration slope mentioned above, and a resultant flow 
spectrum is derived, as shown in Figure . These form observational 
upper bounds on the magnitude of turbulent flows in the convection 
zone at the scales to which the measurements are sensitive. 

It is seen that constraints in Figure become poorer with greater 
imaging depth. This may be attributed to diffraction, which limits 
seismic spatial resolution to approximately a wavelength. In turn, 
the acoustic wavelength, proportional to sound speed, increases with 
depth. Since density also grows rapidly with depth, the velocity re- 
quired to transport the heat flux of a solar luminosity decreases, a 
prediction echoed by all theories of solar convection. Thus we may 
reasonably conclude that the r / Rq — 0.96 curve is also the up- 
per bound for convective velocities at deeper layers in the convective 
zone (although the constraint at r / Rq — 0.92 curve is weaker due 
to a coarser diffraction limit). Less restrictive constraints obtained at 
depths r /Rq — 0.79, 0.86 (whose quality is made worse by the poor 
signal-to-noise ratio) are not displayed here. 

Discussion 

Convective transport. The spectral distribution of power due to an 
ensemble of convective structures, of spatial sizes small or large or 
both, will be broad. For example, it has been argued ([8]) that pho- 
tospheric convection comprises only granules and supergranules, and 
that the power spectrum of an ensemble of these structures would 
extend from the lowest to highest £. In other words, if granulation- 
related flow velocities were to be altered, the entire power spectrum 
would be affected. Thus the large scales which we image here (i.e., 
power for low £), contain contributions from small and large struc- 
tures alike, and represent, albeit in a complicated and incomplete 
manner, gross features of the transport mechanism. 

Our constraints show that for wavenumbers £ < 60, flow veloc- 
ities associated with solar convection (v/Rq = 0.96) are substan- 
tially smaller than current predictions. Alternately one may interpret 
the constraints as a statement that the temporal coherence of con- 
vective structures is substantially shorter than predicted by current 
theories. Analysis of numerical simulations ([6]) of solar convection 
shows that a dominant fraction (~ 80%) of the heat transport is ef- 
fected by the small scales, However, our observations show that the 
simulated velocities are substantially over-estimated in the wavenum- 
ber band £ < 60, placing in question (based on the preceding ar- 
gument) the entire predicted spectrum of convective flows and the 
conclusions derived thereof. We further state that we lack definitive 
knowledge on the energy-carrying scales in the convection zone. We 
may thus ask: how would this paradigm of turbulence affect extant 
theories of dynamo action? 

For example, consider the scenario discussed by [25], who en- 
visaged very weak upflows, which, seeded at the base of the convec- 
tion zone, grow to ever larger scales due to the decreasing density as 
they buoyantly rise. These flows are in mass balance with cool inter- 
granular plumes which, formed at the photosphere, are squeezed ever 
more so as they plunge into the interior. Such a mechanism presup- 
poses that these descending plumes fall nearly ballistically through 
the convection zone, almost as if a cold sleet, amid warm upwardly 
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diffusing plasma. In this schema, individual structures associated 
with the transport process would elude detection because the upflows 
would be too weak and the downflows of too small a structural size 
(M. Schussler, private communication, 201 1). When viewed in terms 
of spherical harmonics, the associated velocities at large scales (i.e., 
low £), which contain contributions from both upflows and descend- 
ing plumes, would also be small. Whatever mechanism may prevail, 
the stability of descending plumes at high Rayleigh and Reynolds 
numbers and very low Prandtl number is likely to play a central role 
([25,3]). 



Differential Rotation and Meridional Circulation. Differential 
rotation, a large-scale feature (£ ~ 2), is one individual global flow 
system and easily detected in our travel-time maps. Differential rota- 
tion is the only feature we "detect" within this wavenumber band. In 
other words, upon subtracting this £ — 2 feature from the travel-time 
maps, the variance of the remnant falls roughly as T _1 , where T is 
the temporal averaging length, suggesting the non-existence of other 
structures at these scales. Consequently, we may assert that we do 
not see evidence for a "classical" inverse cascade that results in the 
production of a smooth distribution of scales. 

Current models of solar dynamo action posit that differential ro- 
tation drives the process of converting poloidal to toroidal flux. This 
would result in a continuous loss of energy from the differentially 
rotating convective envelope and Reynolds' stresses have long been 
thought of as a means to replenish and sustain the angular velocity 
gradient. The low Rossby numbers in our observations indicate that 
turbulence is geostrophically arranged over wavenumbers £ < 60 
at the depth r /Rq — 0.96, further implying very weak Reynolds 
stresses. Because flow velocities are likely to become weaker with 
depth in the convection zone, the Rossby numbers will decrease cor- 
respondingly. At wavenumbers of £ ~ 2, the thermal wind balance 
equation describing geostrophic turbulence likely holds extremely 
well within most of the convection zone: 



on C 

oz r z sin 



where VLq is the mean solar rotation rate, VL is the differential rota- 
tion, z is the axis of rotation, is the latitude, C is a constant, S is 
the azimuthally and temporally averaged entropy gradient. Differen- 
tial rotation around £ ~ 2 is helioseismically well constrained, i.e., 
the left side of equation (2) is accurately known (e.g., [26]). The iso- 
rotation contours are not co-aligned with the axis of rotation, yielding 
a non-zero left side of equation (2). Taylor-Proudman balance is bro- 
ken and we may reasonably infer that the Sun does indeed possess a 
latitudinal entropy gradient, of a suitable form so as to sustain solar 
differential rotation (see e.g., [27, 28]). 

The inferred weakness of Reynolds stresses poses a problem to 
theories of meridional circulation, which rely on the former to effect 
angular momentum transport in order to sustain the latter. Very weak 
turbulent stresses would imply a correspondingly weak meridional 
circulation (e.g., [29]). 
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Fig. 4. Because wavelengths of helioseismic waves may be comparable to or larger 
than convective features through which they propagate, the ray approximation is in- 
accurate and finite-wavelength effects must be accounted for when modeling wave 
propagation in the Sun [20]. In order to derive the 3D finite-frequency sensitivity 
function (kernel) associated with a travel-time measurement [21], we simulate waves 
propagating through a randomly scattered set of 500 east- west-flow 'delta' functions, 
each of which is assigned a random sign so as not to induce a net flow signal [22] (up- 
per panel). We place these flow deltas in a latitudinal band of extent 120° centered 
about the equator, because the quality of observational data degrades outside of this 
region. We perform six simulations, with these deltas placed at a different depth in 
each instance, so as to sample the kernel at these radii. The bottom four panels show 
slices at various radii of the sensitivity function for the measurement which attempts 
to resolve flows at r / Rq = 0.96. Measurement sensitivity is seen to peak at the 
focus depth, a desirable quality, but contains near-surface lobes as well. Note that 
the volume integral of flows in the solar interior with this kernel function gives rise to 
the associated travel-time shift, which explains the units. 



Fig. 5. Observational bounds on flow magnitudes and the associated Rossby num- 
bers. Panels a, b: solid curves with 1-a error bars (standard deviations) show ob- 
servational constraints on lateral flows averaged over m at radial depths, r / Rq = 
0.92, 0.96; dot-dash lines are spectra from ASH convection simulations [6]. Colours 
differentiate between the focus depth of the measurement and coherence times. At 
a depth of r / Rq = 0.96, simulations of convection [6] show a coherence time of 
— 24 hours (panel a) while MLT [16] gives = 96 hours (panel b), the 

latter obtained by dividing the mixing length by the predicted velocity. Both MLT 
and simulations [23, 24] indicate a convective depth coherence over 1.8 pressure scale 
heights, an input to our inversion. At r / Rq = 0.96, MLT predicts a 60 ms _1 , 
t = 61 convective flow and for r / Rq = 0.92, an £ = 33, 45 ms — 1 flow (upon 
applying continuity considerations [23]). Panel c shows upper bounds on Rossby num- 
ber, Ro = U/(2QL), L = 27rr/y/e(£+l), r = 0.92, O.96i? . Interior 
convection appears to be strongly geostrophically balanced (i.e., rotationally domi- 
nated) on these scales. By construction, these measurements are sensitive to lateral 
flows i.e., longitudinal and latitudinal at these specific depths (r /Rq = 0.92, 0.96) 
and consequently, we denote these flow components (longitudinal or latitudinal) by 
scalars. 
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